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Abstract 

Energy spectrum of cosmic rays (CR) exhibits power-like behavior with very char- 
acteristic "knee" structure. We consider a generalized statistical model for the pro- 
duction process of cosmic rays which accounts for such behavior in a natural way 
either by assuming existence of temperature fluctuations in the source of CR, or by 
assuming specific temperature distribution of CR sources. Both possibilities yield 
the so called Tsallis statistics and lead to the power-like distribution. We argue that 
the "knee" structure arises as result of abrupt change of fluctuations in the source 
of CR. Its possible origin is briefly discussed. 
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1 Introduction 



The origin of the characteristic features of the energy spectrum of cosmic rays (CR), 
which has power-like behavior with "knee" structure, remains matter of hot debate 
(for survey of models proposed to explain the origin of CR see []]). It could reflect 
different regimes of diffusive propagation of CR in the Galaxy, but it could also 
be due to some property of acceleration processes within the source of the CR 

* Corresponding author 

Email addresses: wilk@fuw.edu.pl (G. Wilk), wlod@pu.kielce.pl (Z. Wlodarczyk). 



Preprint submitted to Elsevier 



29 October 2008 



itself. In this second crucial question is whether the sources of CR below 

the "knee" can also accelerate particles to much higher energies, so that a single 
population of astrophysical objects can explain the smooth spectrum of cosmic rays, 
as observed over many orders of magnitude in energy. We address this problem using 
a generalized statistical model specially adapted to this end. However, our work will 
concentrate more on the physics of CR than on the generalized statistics, providing 
therefore some physical explanations to ideas presented already in [2] and [3J. In 
particular, we shall argue that the observed "knee" structure of the CR energy 
spectrum arises as result of some abrupt change of fluctuation pattern in the source 
of CR. 



2 Nonextensive statistics and results 



We shall start with some basic information on nonextensive statistical mechanics as 
introduced by Tsallis [I], which have been already successfully applied to a variety of 
complex physical systems, including CR where, among others, the energy spectrum 
of cosmic rays have been analyzed from a nonextensive point of view [2|3] . The idea 
is to maximize the more general entropy measure than the usual Boltzman-Gibbs- 
Shanon (BGS) entropy, the one which depends on a new additional parameter q and 
which leads to generalized version of the statistical mechanics. If we optimize, under 
appropriate constrains, the BGS entropy, 

S = - J dEP(E) In P(E), (1) 
we obtain the equilibrium distribution in the usual form of exponential distribution, 

P(E) = \ exp (-|) . (2) 

This distribution can alternatively be obtained as the solution of simple differential 
equation: 

dP{E) P(E) 

~1e~ ~ ~T~- [6) 

A more general formalism proposed in [I] (and sometimes referred to as nonextensive 
statistical mechanics) is based on the generalized entropy, 



Its maximization under appropriate constrains yields a characteristic power-like dis- 
tribution (which sometimes is also called g-exponential distribution, exp g (. . . )): 
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For g — > 1 one recovers the usual exponential distribution This equilibrium dis- 
tribution can alternatively be obtained by solving the following differential equation 

dP(E) = _PW 
dE T ' U 

Precisely this equation has been used in [2] to describe the flux <&(E) of cosmic rays. 
However, the values of the temperature obtained there seem to be uncomfortably 
high. 

On the other hand there is growing evidence that nonextensive formalism applies 
most often to nonequilibrium systems with a stationary state that posses strong 
fluctuations of the inverse temperature parameter (3 = 1/T [5116] . In fact, fluctuating 
(3 according to gamma distribution with variance Var((3) results in a power like 
distribution (JSJ) with nonextensivity parameter q being given by the strength of 
these fluctuations, 

distribution (j2J) is just its limiting case when q — > 1. This observation was used 
in [3] to describe the flux <&(E). Although the results were reasonably good the 
estimated temperature T ~ 170 MeV (comparable with the so called Hagedorn 
temperature known from description of hadronization processes) seems to be, again, 
overestimated. This was because author insists on description of the whole range of 
energy spectrum including its very low energy part which, in our opinion, is governed 
mainly by the geomagnetic cut-off and should be considered separately. 



2.1 Energy spectrum 



For relativistic particles (where the rest mass m can be neglected) the energy E ~ p 
and the density of states of an ideal gas in three dimensions is given by u(E) oc E 2 . 
The flux <&(E) can be then obtained straightforwardly from P(E) and reads 

= N E 2 P(E) } (8) 

where No is normalization factor. For E » T we have power spectrum $(£') oc E^ 1 
with the slope parameter 7 which in terms of parameter q introduced above is 
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7 = (3 — 2q) /(q — 1). In the case of CR this spectrum has the shape of broken power 
law E^ 1 which changes pattern in the region named as "knee" with the slope 71 ~ 2.7 
at energies below ~ 10 15 eV and 72 — 3.1 at energies between the "knee" and the 
highest measurable energies E ~ 10 18 eV. In the language of the nonextensivity 
parameters it would mean that qi = 1.213 before and q2 = 1.196 after the "knee", 
i.e., one can argue then that at the "knee" one witnesses the change of fluctuation 
pattern. 

2.2 Temperature fluctuations 

As mentioned above the special role in converting exponential distribution to its q- 
exponential counterpart play fluctuations of the scale parameter provided in the form 
of gamma function. There are at least two scenarios leading to gamma distribution 
in ft mentioned above (here ^T 1 = ft (q — 1) and v~ l = q — 1): 



(a) temperature distribution of sources; 

(6) temperature fluctuations in small parts of a source. 

In what concerns the first possibility notice that gamma distribution ([9]) is the most 
probable outcome of the maximalization of Shannon information entropy (JTJ) under 
constraints that / f(ft)dft = 1, (ft) = fto and (because distribution we are looking 
for is one sided, i.e., defined only for ft > 0) that (ln(/5)) = ln(z//3 ). However, this 
possibility is rather unlikely because in this case one expects that there is some cut 
temperature T cut such that T < T cut , which would result in very characteristic rapid 
break in the CR energy spectrum, not observed in experiment. 

To illustrate the second possibility let us suppose that one has thermodynamic 
system, different small parts of which have locally different temperatures, i.e., its 
temperature understand in the usual way fluctuates. Let £(t) describes the stochas- 
tic changes of temperature in time and let it be defined by the white Gaussian noise 
((£(£)) = and (£(t)£(t + At)) = 2D5(At)). The inevitable exchange of heat which 
takes place between the selected regions of our system leads ultimately to the equi- 
libration of temperature. As we have advocated in [5], the corresponding process 
of the heat conductance leads eventually to the gamma distribution f(ft) (jSJ) with 
variance ([7]) related to the specific heat capacity Cy of the material composing this 
system by 




(9) 




(10) 
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The change of fluctuation pattern in the "knee" region mentioned before would 
therefore correspond in this case to abrupt change in the heat capacity of the order 
of C 2 /d = 1.09. 
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Fig. 1. Example of the possible modification of gamma distribution needed to describe 
fluctuation of temperature which could lead to description of the "knee region" of the CR 
spectra. At T = 10 9 MeV parameter q changes from q = 1.214 (it would result in dotted 
line) to q = 1.2 (in this case curve changes slope and continues as full line). Such break 
represents some abrupt change in the fluctuation patter, which we attempt to explain here. 

2.3 Heat capacity 



Can one expect something of this kind to happen in the astrophysical environment 
of the CR? In what follows we shall argue that, indeed, one can. Let us first notice 
that subject of temperature fluctuations in astrophysics is much-discussed problem 
nowadays. Its effect on the temperatures empirically derived from the spectroscopic 
observation was first investigated in [7] whereas in [8|9|10|llj it was shown that 
temperature fluctuations in photoionized nebulea have great importance to all abun- 
dance determinations in such objects. It means that discussion of the heat capacity 
or, equivalently, the behavior of parameter q defining the energy spectrum, is fully 
justified. 

Let us concentrate therefore on the problem of heat capacity of astrophysical objects, 
in particular in neutron stars. In such objects one observes the following feature. The 
total specific heat of their crust, C, is the sum of contributions from the relativistic 
degenerate electrons, from the ions and from degenerate neutrons. In the tempera- 
ture that can be reached in the crust of an acreating neutron star (which is of the 
order of T ~ 5 • 10 8 K and is below the Debaye temperature Tp ~ 5 ■ 10 9 K) we have 
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Ci on < C e < C n . When the temperature drops below the critical value T = T C the 
neutrons become superfluid and their heat capacity C%f increases [T2][T3] . 



c„ 



9 1 V^C 
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+ 3.68 



T 



(11) 



At T ~ 0.7Tc we have Cff ~ l.lC n what corresponds to the changes of spectral 
index by A7 ~ 0.5. To summarize: one witnesses here the abrupt change in the heat 
capacity at some temperature. 
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Fig. 2. CR energy spectra fitted by simple Tsallis distribution (left panel) and distribution 
using specially adopted fluctuation function example of which is presented in Fig. []] (right 
panel). Notice the total inadequacy of simple exponential (Boltzman, denoted by BG here) 
distribution in description of these spectra. 

Suppose now that we take seriously the conjecture expressed by eq. (TlOT) that Cy is 
directly connected with the parameter q. It means then that, the usual fluctuation 
pattern given by gamma distribution should be replaced by it slightly modified 
version shown in Fig. [TJ, which is characterized by two nonextensivity parameters, q± 
and <72- This change is assumed to be abrupt and taking place at some temperature 
T cut and it differs our proposition from what was proposed in [2], which from our 
point of view, corresponds to some form of composition of two (suitably normalized) 
such gamma functions with different q each. Following our proposition one obtains 
the following flux of CR: 



$(E) = N E 2 - [P gi (E) - a 1 {E)P gi {E) + a 2 (E)P g2 (E)} 

where P q (E) is given by eq. (JSD and 

1 1 - (1 - qi ) E/T 



(12) 



qi-1 (q { -l)T cut /T 



/r 
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Our results are presented in Fig. [2J The "knee" region is reproduced very well, 
however, the price to be paid is the need of suitable choice of energy at which 
fluctuation pattern changes (cf. Fig. [TjQ. 



2.4 Acceleration 



As one can see in Fig. [2] we have obtained agreement with data for the whole range of 
CR energy spectrum but the price paid for this is, again, apparently too high value 
of the temperature parameter used, T = 100 MeV. We have prone therefore to the 
same kind of criticism as we have applied to previous attempts in this field [2|3] . The 
possible way out of this dilemma is to argue that the break in the original spectrum 
and connected with it phase transition occurs actually at much slower energies and 
that resultant spectrum is then accelerated to the observed energies - for example 
by the magneto-hydrodynamical turbulence and/or shock discontinuities (i.e., by 
the so called DSA mechanism, cf. [19]). The simplest version of this mechanism, 
as discussed in [20], implies that distribution function after the shock, f a fter{p), is 



related to the original distribution before the shock, fb e fore{p) in the following way: 

fafterip) = ~ b f <W M f before (p'), (13) 



where b = 3r/(r — 1). Here p denotes the particle momentum, r = P2/P1 describes 
the compression of densities across the shock and p m in denotes the minimal value of 
momenta. DSA mechanism transforms a 5 (p — po) spectrum of relativistic particles 
in a power-like spectrum of the type f(p) oc p~ b . For example, if one has initial 
spectrum of the form oc p~ c which encounters a shock with strength given by b then 
eq. f|T3|) shows that: 



1 Two remarks are in order at this point. First, the nucleon superfluidity was predicted 
already in [T3] and today pulsar glitches provide strong observational support for this hy- 
pothesis |15| . Nucleon superfluidity arises from the formation of Cooper pairs od fermions 
(actually in [16] also quark superfluidity from cooling neutron stars were investigated). 
Continuous formation and breaking of the Cooper pairs takes place slightly below T = Tq 
(critical temperature Tq is in the order 10 9 — 10 10 K). The other, neutron stars are born 
extremely hot in supernova explosions, with interior temperatures around T ~ 10 12 K. 
Already within a day, the temperature in the cental region of the neutron star will drop 
down to oc 10 9 — 10 10 K and will reach 10 7 K in about 100 years [T7]. The first measure- 
ments of the temperature of a neutron star interior (core temperature of the Vela pulsar 
is T ~ 10 8 K, while the core temperature of PSR B0659+14 and Geminga exceeds 2 • 10 8 
K) allow to determine the critical temperature Tc ~ 7.5 • 10 9 K [18] . 
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in the case of b > c (which corresponds to the initial spectrum being softer than 
it would result from a 5-function injected into the shock) one has f after oc P c , 
i.e., the acceleration does not change the shape of the spectrum; 
in the case of b < c (i.e., for the steep initial spectrum) one has f after p~ b , which 
coincides with result of injection of a ^-function into the shock. 



It can be then shown from eq. ([13]) that when an ensemble of shocks is encountered, 
the shape of the spectrum should be given by the strongest shock [20] . 



That such scenario is a priori plausible in the case considered here can be seen by 
the following argumentation. If the energy growth is given by 



dt 

than from master equation 



= a + bE, (14) 



9 -T - < i5 > 

one gets the evolution equation 

dP(E) dt P(E) 

-^E- = - cPiE) dE = -T(E) { ) 

in the form of eq.([3]) with energy dependent temperature parameter: 

T( E) = = ^„ + (,-!)£ 

c q 

where we have used: c = T _1 , a = q~ l and b = (1 — g _1 )T _1 . Energy dependent 
temperature parameter T(E) in this form immediately leads to the energy distribu- 
tion given by eq.flSD- Notice that q^ 1 plays here the role of the weight with which 
we select the constant (thermal) and linear (accelerating) terms in the equation de- 
scribing the growth of energy. It means therefore that Tsallis distribution preserves 
its structure when subjected to the aforementioned acceleration scheme. The prob- 
lem, which remains to be solved is whether the broken spectrum with the "knee" 
structure can be suitably transformed in the same way. In particular, whether the 
"knee" structure is preserved and how its position before the acceleration process is 
related to that actually observed. 

To summarize this part: stochastic mechanisms of acceleration of CR particles (like 
acceleration on the fronts of shock waves or Fermi acceleration in turbulent plasmas, 
both analogous in some sense to Brownian motion) do not change the shape of the 
production spectra, but, unfortunately, they are not particularly effective, i.e., they 
do not lead to large increase of energy [21]. The increase of energy per one collision 
is of the order AE/E ~ u 2 /c 2 , what for the plasma velocity u ~ 10 6 — 10 7 cm/s 
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gives AE/E w and leads to the mean relative increase of energy during the 
time life of Galaxy (t « 3 • 10 17 s) only by factor (5) « 3. Fluctuation on the steep 
spectrum of accelerated particles result in additional increase of energy. Because 
of the multiplicative character of acceleration we have log-normal distribution of 
variable 5, P(5)d\n5 = (v 7 ^)" exp [- (In 5 - ln(<5)) 2 /2 a 2 \ din 5, what results in 
shift of the spectrum of source on the energy scale by (£)(t- 1 Vt ex p _ l) 2 a 2 / (27)], 
where a 2 is variation of the distribution P(8). For 7 w (5) « a 2 « 3 we can obtain 
only order of magnitude shift of the energy spectrum. 



3 Summary and conclusions 



Let us summaries arguments presented above. 

• The spectrum of cosmic rays has the shape of broken power law E' 1 , with the 
slope 71 ~ 2.7 at energies below ~ 10 15 eV and 72 ~ 3.1 at energies between the 
knee and E ~ 10 18 eV. This slopes correspond to the nonextensivity parameters 
(taking into account that q = (3 + 7)/(2 + 7)) qi = 1.213 and q 2 = 1.196, 
respectively, i.e., to the change of heat capacity of the order of C 2 /Ci = 1.09. 

• Nonextensive statistics successfully describes the smooth power-law spectrum and 
traces its origin back to fluctuations of temperature, f(/3), being given by gamma 
distribution. 

• Out of two possible scenarios leading to such distribution of inverse temperature 
we prefer the temperature fluctuation in the source rather than the temperature 
distribution of sources. The point is that in the second case some cut temperature 
T cut is expected which would result in the rapid break in the energy spectrum 
and which is not observed. The temperature T (not essential in the high energy 
region, E » T) seems to be of the order of MeV, i.e. of the order of the interior 
stars temperature (if stars are born extremely hot in supernova explosions, with 
interior temperatures around T ~ 100 MeV, already within a day the temperature 
in the cental region of the star will have dropped down to ~ 0.1 — 1 MeV and 
reach the 1 keV in about 100 years). The critical temperature (corresponding to 
the nucleon superfluidity) is Tq ~ 0.1 — 1 MeV. It means then that the origin of 
changes of the nonextensivity parameters at temperature T cut ~ 10 15 eV ~ 10 19 
K is still open question. 

• The nonextensive formalism leads to production (injection) spectrum and the 
acceleration processes (dE/dt ~ E, which does not change the shape of power 
spectrum) and allows the shift of this spectrum to high energies. The question of 
how to connect "knee" position before and after such acceleration process remains, 
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however, still open. 



Let us close with some very intriguing observation. Independently of the discussion 
presented above one can notice that the measured CR energy spectrum can be 
converted using arguments presented above (eq. fflUj) connecting Tsallis parameter q 
with the heat capacity) into energy dependence of the heat capacity C. The result is 
shown in Fig. [3j As one can see C acts here as a kind of magnifying glass converting 
all subtle structures of <&(E) into much more pronounced and structured bump. 
Its importance would parallel long-standing discussion of the origin of the knee-like 
structure of the CR energy spectrum, but exposed in much more dramatic and visible 
way. At the moment we can only offer two examples of the possible explanation of 
this feature. The first is that this effect is due to the change of the effective number 
of degrees of freedom in the incoming projectiles with energy. Assuming that most 
of CR consist of protons, which are build from three quarks, one could speculate 
that each bumps correspond to excitations from single proton to proton plus one 
quark and two quark structure, after exciting all three quarks one comes back to 
the original situation (much in the spirit of the changes observed when ice becomes 
water and this becomes steam - there also corresponding C"s show characteristic 
jump [22]). The other, perhaps more realistic explanation of Fig. El is to connect it 



5.3 
5.2 
5.1 



4.8 

4.7 1 — 1 

1E9 1E101E11 1E121E131E141E151E161E171E181E191E20 

E[eV] 

Fig. 3. Pattern of the energy dependence of the heat capacity obtained from the CR energy 
spectrum in the vicinity of the "knee" structure using nonextensive statistical approach. 

with the behavior of heat capacity in Fermi liquids. Following [12] the proton heat 
capacity C ~ m*/m, i.e., it is proportional to the ratio of the effective mass of the 
proton in the neutron fluid to the mass of the free proton. In the case of a mixture of 
Fermi liquids the proton effective mass m* is affected by interactions with neutrons 
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and other protons and is given by 




^ = 1+ l D *> /r+(^l /r , as) 

where D p denotes the density of quasiparticle states at Fermi surface given by wave 
vectors k Fn and k Fp for, respectively, neutrons and protons whereas f pp and /f n are 
Landau parameters [23]. Fig.[3]can be then interpreted as showing changes of C with 
energy in the Fermi liquid. We start with the superfluid liquid with C\ = 4.7 (here 
m* represents effective mass for pp and pn interactions), when energy increases we 
stop to see nuclear interactions and C 2 = 5.1 (with m* representing pp interactions 
only), finally, for large T, one has Fermi gas with C3 = 5.3 and, still further, the 
usual Fermi liquid^- Notice that 

l n f pp C 2 — C 3 , 1 ^ ( k Fn \ rVn C\ — C 2 



what results in the following relation between Landau parameters, 



2 /r d-c 3 



^f p J fi C 2 — C 3 



2. (20) 



In the case of a one-component Fermi liquid we have the well known identity, 
m* jm = l + .Ff p /3, where Ff p = D p fl P . From (1201) we can see that in two-component 
Fermi liquid the quantity 1 — m*/m is 3 times bigger (this is because parameter 
fi which determines interaction between quasiparticles is negative what results in 
smaller effective mass). From properties of excited states in nuclear matter (Pb and 
neighbor nuclei [25]) F" p = -0.5 ± 0.25. If < F± p < Ff p and taking (after 
[26J) -F™ n — Fi P = —0.2, we can estimate that for neutron-star matter one has 
m*/m = 1 + Ff p ~ 1 - 0.4 ± 0.3 = 0.6 ± 0.3. 

We end by saying that the above are just plausible examples, not fully convincing 
explanation(s). It means then that this problem deserves further scrutiny to be 
performed elsewhere. 



2 It is worth to remember that fluctuations of temperature we are talking about in this 
work refer to fluctuations in small region V. For Fermi liquid the heat capacity expressed 
in units of Boltzmann constant ks (i-e., for kg = 1) is of the order C ~ 3 • 10 35 cnr 3 |24J. 
Therefore, taking values of C estimated from the slope of the primary CR spectra (cf. Fig. 
[3|) one gets that size of the region of fluctuations is V ~ 10 4 fm 3 . 
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